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Abstract 

An explanation is given for previous numerical results which suggest a certain bifurcation of 
'vector solitons' from scalar (single-component) solitary waves in coupled nonlinear Schrodinger 
(NLS) systems. The bifurcation in question is nonlocal in the sense that the vector soliton does 
not have a small- amplitude component, but instead approaches a solitary wave of one component 
with two infinitely far-separated waves in the other component. Yet, it is argued that this highly 
nonlocal event can be predicted from a purely local analysis of the central solitary wave alone. 
Specifically the linearisation around the central wave should contain asymptotics which grow at 
precisely the speed of the other-component solitary waves on the two wings. This approximate 
argument is supported by both a detailed analysis based on matched asymptotic expansions, 
and numerical experiments on two example systems. The first is the usual coupled NLS system 
involving an arbitrary ratio between the self-phase and cross-phase modulation terms, and the 
second is a coupled NLS system with saturable nonlinearity that has recently been demonstrated 
to support stable multi-peaked solitary waves. The asymptotic analysis further reveals that 
when the curves which define the proposed criterion for scalar nonlocal bifurcations intersect 
with boundaries of certain local bifurcations, the nonlocal bifurcation could turn from scalar to 
non-scalar at the intersection. This phenomenon is observed in the first example. Lastly, we 
have also selectively tested the linear stability of several solitary waves just born out of scalar 
nonlocal bifurcations. We found that they are linearly unstable. However, they can lead to 
stable solitary waves through parameter continuation. 
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1 Introduction 

Solitary waves play an important role in the solution dynamics of nonlinear evolution equations. If 
the solitary waves are stable, they often emerge as final states in an initial-value problem. Even if 
these waves are unstable, the mere existence of such waves has important implications for solution 
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evolutions. In recent years, it has been discovered that complicated solitary waves could bifurcate 
from simple solitary waves. 

Local bifurcations of 'wave and daughther waves' ('vector solitons') from single-component waves 
('scalar solitons') have been studied in various forms of coupled nonlinear Schrodinger (NLS) sys- 
tems [4, 5, 6, 10, 14, 17]. The bifurcation is local in that bifurcated solitary wave is infinitesimally 
close to the original solitary wave (as a graph) at the point of bifurcation. The condition for such 
a local bifurcation to occur is based on the linearisation around the single-component pulse having 
a solution with purely decaying asymptotics at infinity (see Sees. 2 and 3 below). 

Nonlocal bifurcations are where the bifurcated solitary wave is not infinitcssimally close to the 
original wave at the point of bifurcation. These have been reported numerically in the usual 
coupled NLS equations (the generalised Manakov system, involving an arbitrary ratio between the 
self-phase and cross-phase studied in Sec. 3 below, henceforth referred to as the CNLS system) 
[5, 14, 15], and in second- and third-harmonic generation systems [6, 18, 19]. If the bifurcated wave 
looks like several vector solitons glued together, that nonlocal bifurcation has also been treated 
analytically by an asymptotic tail-matching method [15] (see Sec. 4 for application of this method 
in the context of this paper). General mechanisms have also been identified that lead to such 
vector non-local bifurcations in classes of coupled NLS systems, such as eigenvalue degeneracy or 
the existence of a local bifurcation [15, 18, 20]. The linear stability of these multiple-pulsed vector 
solitons for the CNLS system has been studied in [16], and it has been shown that such states are 
linearly unstable. However, if the bifurcated wave is glued together by scalar (i.e., single-component) 
solitary waves, no analysis has been performed to our best knowledge. In this paper, we focus on 
this scalar nonlocal bifurcation. We will show that it is closely related to local bifurcations, and 
can be treated on the same footing. 

The criterion we propose for a scalar nonlocal bifurcation (in Sec. 2) is that the solution of the 
linearized equation around the central single-component pulse should have only purely growing 
asymptotics instead of purely decaying asymptotics at infinity. In the following Sections, we test 
this criterion against two example systems. The first, in Sections 3, is the CNLS system. The 
second, in Section 5, is a coupled NLS system with saturable nonlinearity. In both cases, good 
agreement is obtained between our bifurcation condition and the numerics. Near the bifurcation 
point, we have also developed a detailed asymptotic analysis based on the above-mentioned tail- 
matching method, which is performed for the CNLS system in Sec. 4. This analysis produces an 
explicit formula for the spacing between scalar solitons being pieced together, and this formula 
agrees well with the numerics. An extension of that analytical theory to more general coupled NLS 
systems, such as our second example studied in Sec. 5 is straightforward, and will be omitted in 
this paper. Lastly, we discover in the course of our asymptotic analysis that if the curve which 
defines our proposed criterion for scalar nonlocal bifurcations intersects with boundaries of certain 
local bifurcations, the nonlocal bifurcation may turn from scalar to non-scalar at the intersection. 
This phenomenon indeed occurs in the CNLS system. We have also studied the linear stability of 
solitary waves just born out of scalar nonlocal bifurcations. The results suggest that these waves 
are always linearly unstable. However, they can lead to stable solitary waves through parameter 
continuation in the second model. 



2 



2 A geometric argument 



Consider a general system of ordinary differential equations (ODEs) of the form 

Uxx-u + fi{u,v)u = 0, (2.1) 

Va;x - u;'^v + /2(n, v)v = 0, (2.2) 

where < a; < 1. It is assumed that /i and /2 are smooth nonlinear functions of their arguments 

which vanish as {u, v) (0, 0) and may well depend on other system parameters. Moreover they 
are such that the problem with ?; = or with v = 0, which are both invariant subspaces within the 
four-dimensional phase space of the ODE system (2.1), (2. 2), contain even homoclinic solutions 

{u{x),v{x)) = {ufi{x),0), Ufi{—x) = Ufi{x), u — >■ as a; — > ±oo, (2-3) 

{u{x),v{x)) = {0,Vh{x)), Vhi-x) = Vh{x), v^O as x ^ ±oo. (2.4) 

In what follows we shall use the terms 'pulse' 'homoclinic to the origin', 'solitary wave' and 'soliton' 
entirely synonymously. We shall refer to the invariant subspace homoclinic solutions (2.3) and (2.4) 
as scalar u- and u-pulses respectively. These scalar solitons are contrasted with vector solitons which 
are homoclinic solutions that have non-zero u and v components. Note that all homoclinic solutions 
to the origin are generic, that is they persist under parameter perturbation, since the system is 
both reversible and Hamiltonian [3]. 

Now consider the linearisation of (2.1), (2. 2) around the n-pulse 

Uxx = U+ ^-^fi[uh{x),{)\uh{x) + fi[uh{x),{)\^u, (2.5) 

Vx:,=Uj\ + f2[uh{x),{))]v. (2.6) 

Note that the linear equations decouple. So let us look at the specific class of solutions to this 
linear problem which have n = 0. Now we have to simply solve the second equation (2.6). The 
general asymptotics of such solutions satisfy 

V cfe-^'' + c^e'^^ + o (e^^'l ''!) as \x\ ±oo (2.7) 

for some constants and c^. Note that we arc able to write the o(-)-term by the assumption that 
< a; < 1 so that the asymptotics of /2(u(x),0) decay more rapidly than exp(— ci;|x|). 

Consider even solutions (2.7) of (2.6), hence = c~ := Ci, i = 1,2. This defines a unique solution 
(up to scale) for the linear initial-value problem (2.6) for x > 0. Now suppose that at a particular 
value of u>, this solution had a particular tail asymptotics (2.7) with C2 = (see Fig. 1(b)). Then 
the solution to the linear problem would be localized. Going back to the fully nonlinear problem, by 
standard bifurcation theory results, we have satisfied the necessary condition for the /oca/ bifurcation 
of a wave and daughter wave consisting of the mother tt-pulse and a small-amplitude i;-component. 
The vanishing of C2 is a codimension-one condition, and hence local bifurcations will lie on lines in 
a parameter plane. As already mentioned in the Introduction, the existence of such bifurcations in 
coupled NLS systems have been established by a number of authors. 

Suppose instead that we find a solution at some parameter value that satisfies ci = 0. Note 
that this describes a pure exponentially growing solution to leading order as x ^ ±00. See 
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Fig. 1(a). However, consider the nonlinear implications of this within the four-dimensional phase 
space {(u,u' ,v,v') : u = u' = 0} of the ODE (2.1), (2. 2). We have found an initial condition for 
an even solution that is an infinitesimal perturbation of Uh{0), and which is attracted as x ^ oo 
towards the invariant plane V := {{u,u' ,v,v') : u = u' = 0}. Moreover, this rate of attraction 
is faster (~ exp(x)) than the exponential contraction or expansion with that plane (~ exp(±c<;a:;)) 
near the origin. Hence the condition that < a; < 1 ensures that the invariant plane V is normally 
hyperbolic. Using standard results for normally hyperbolic manifolds, the eventual behaviour of 
this perturbed trajectory is governed by its behaviour on V . The fact ci = implies that the 
trajectory is attracted onto the local unstable manifold within V . But the local unstable manifold 
on V is precisely the piece of trajectory that forms the u-pulse solution Vh. Also since we are talking 
about an infinitesimal perturbation to the underlying pulse Uh, the time taken to be attracted to 
Vh in this way is arbitrarily long. Hence we have the scalar nonlocal bifurcation of two Vh pulses 
'at ±oo' as depicted in Fig. 1(a). 

The above is only a plausible argument, but it is highly appealing from an intuitive point of view. 
In Sec. 4, we will develop an in-depth asymptotic theory for this scalar nonlocal bifurcation in the 
CNLS equations. The results of that theory fully support the above intuitive argument. 

Before proceeding to the examples, let us make a few short remarks here. 

1. First, note the need to assume that < u; < 1. This was required in order to make the 
u = Q invariant plane normally hyperbolic, or equivalently to assume that the asymptotic 
attraction onto this plane was o{e~'^^). If this condition is violated, then there is no sense in 
which the perturbed trajectory converges only to the unstable manifold within V and hence 
the argument fails. 

2. Second, it is interesting to note that our analysis suggest that local and scalar nonlocal 
bifurcations can be treated on an equal footing. One requires that ci vanishes, the other that 
C2 vanishes. It is perfectly possible to imagine a scenario where, as a parameter is varied, the 
even solution v{x) to the linear problem (2.6) generates extra internal oscillations. If it does 
so in a smooth way, then it is clear that it must pass repeatedly through successive zeros of 
ci and C2. Hence one would find each scalar nonlocal bifurcation sandwiched between two 
successive local bifurcations. 

3. In fact, in the above it is not quite enough to assume that ci vanishes for a scalar nonlocal 
bifurcation. We must have that C2 has the correct sign to attach to the component of the 
unstable manifold of the nonlinear equation that contains the pure u-pulse. In all the examples 
below, the pure-?; equation is odd and hence both and —Vh are solutions. Thus either sign 
of C2 will lead to a scalar nonlocal bifurcation. 

4. Finally, in systems with odd nonlinearity, the above arguments can be repeated to find anti- 
symmetric scalar nonlocal bifurcations where the two daughter waves at infinity are and 

We now turn to two examples to test the validity of this approximate reasoning. 
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3 Example 1: the CNLS equations 



The usual coupled nonlinear Schrodinger (CNLS) equations may be written in dimensionless form 
as 

iUt + U^^ + m'' + (3\V\'')U = 0, (3.1) 
iVt + V^^ + {\V\'' + l3\U\'')V = 0. (3.2) 

They have been used to describe the interaction between wave packets in dispersive conservative 
media, and also the interaction between orthogonally polarized components in nonlinear optical 
fibres; see [14, 5] and references therein. Looking for steady solutions of the form 

[/ = e^'^^Va;), V = e'^'^^vix), 

and performing scaling so that = 1 and 0^2 = we arrive at the following set of ODEs 

Uxx — u + {u^ + (5v^)u = 0, (3.3) 
- u?v + {v'^ + (iu^)v = 0. (3.4) 

Here /? is a real and positive cross-phase-modulational coefficient, and lo (> 0) is a propagation 
constant parameter. It is noted that if [u{x;uj),v{x;(jj] is a solution, then another solution at 
propagation constant ^ can be obtained via the transformation [14] 

[-v{-; w), -u{-; w)]. (3.5) 

UJ U> LO LO 

Thus, in this paper, we restrict < a; < 1. 

We note that when (3=1, the partial differential equations (PDEs) (3.1), (3.2) are called the 
Manakov system which is integrable [8]. In this case, all solitary waves of the ODE system (3.3) 
and (3.4) have closed-form analytical expressions [1]. When /? = 0, the PDEs arc two copies of the 
single NLS equation which is also integrable [21]. The solitary waves for /3 = are simply sech 
pulses. When /3 7^ or 1, the structure of solitary waves in this system is much more complicated. 
This structure was partially unraveled in [4, 2, 1, 5, 14, 15]. It is known that local bifurcations 
occur along curves in the (/3, a;)-plane that are given by closed form expressions (sec below). These 
local bifurcations are where wave and daughter-wave structures are born. In other words, at local 
bifurcations, a small and localized v component develops from a pure -u-pulse. It was also observed 
numerically that scalar nonlocal bifurcations such as described in this paper occur. That is, passing 
through the bifurcation is a single component pulse, for which v is zero everywhere, and u is given 
by a sech function. Bifurcating from this is a solution for which the u component remains about the 
same, but the v component suddenly develops two pulses which are far separated from the central 
u pulse. These v pulses can be symmetrically or antisymmetrically distributed. Their sizes jump 
from zero to a certain finite size across the bifurcation. In this section, we analytically determine 
the boundaries of these scalar nonlocal bifurcations through the criterion developed in Sec. 2 and 
compare them with direct numerical results. 

3.1 Local and scalar nonlocal bifurcations 

First, we recall the results for local bifurcations in this system [4, 17, 14], assuming that a small 
?;-component bifurcates from a pure -u-pulse. Thus, at the bifurcation point, the v component is 



5 



infinitesimally small. Thus the ti-component is simply governed by the equation 

Uxx — u + = (3-6) 

whose homoclinic solution is 

u{x) = V2sech.x. (3.7) 

According to standard results, a necessary condition of a local bifurcation of a homoclinic solution 
with a small-amplitude v component from the u-pulse (3.7) is that there is a non-trivial localised 
solution to the linearised problem of the t;-component. This takes the form of a linear Schrodinger 
equation 

Vxx - + 2/3sech2a;u = 0, (3.8) 

and for local bifurcation we require 

— > as \x\^ iboo. 

This equation can be solved exactly [7] as we now explain. With the variable transformation 

V = sech^x '0, ^ = sinh^ x, (3.9) 

where 



s = , (3.10) 



the Schrodinger equation (3.8) becomes 



e(l + 04>xx + [(1 - s)i + ]^](t>x + \{s''- UJ^)<t> = 0, (3.11) 



which is a hyper-geometric equation. Its even and odd solutions are 

4>i = - \s, - \s, \, -0, (3.12) 

</'2 = VlF{]^u - Is + 1, -i., - Is + 1, 1 -^), (3.13) 

where F is the hyper-geometric function. In order for the solution vi = sech^x^i to decay to zero 
as X goes to infinity, one must have 

11 , , 

2^~2^ = -'^b (3-14) 

where n\ is a non-negative integer, and n\ < ^s. Then the v\ solution is 

k=Q K2)k'^- 

which decays to zero as ^"'^/^ (i.e., e"'^^). In the solution (3.15), {a)k is defined as 

,.\a{a + l){a + 2)...{a + k-l), k>0 

{a)k = < T (3.16) 



In order for solution V2 = sech^x4>2 to decay to zero as x goes to infinity, one must have 

111 ^ ^ 

-W - -S - = -712, (3.17) 
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where n2 is a non-negative integer, and n2 < ^(s — 1). Then the V2 solution is 

V2 = sech a; smh x } ^ 3 , (3.18) 

jk=o (2)^^' 

which also decays to zero as ^""'^/^ (i.e., e^^''^). 

When conditions (3.14) and (3.17) are combined, we find that the boundaries for local bifurcations 
are 

u = Lo^^{l3) = s-n, (3.19) 

where s is given by (3.10), n is a non-negative integer, and n < s. The first boundary (n = 0) 
exists for any /? > 0; the second boundary (n = 1) exists only for /? > 1; the third boundary 
(n = 2) exists only for /? > 3: etc. The first three boundaries ajg ^2 plotted in Fig. 2 as dashed 
lines for illustration. Note from the above construction that even n corresponds to the existence of 
symmetric bifurcating waves (even in both u and v) whereas odd n corresponds to anti-symmetric 
bifurcation (even in u, odd in v). 

Now, how can we define the scalar nonlocal bifurcations? It can be noted that on the above local 
bifurcation boundaries, the appropriate solution vi or V2 has the following asymptotic behavior: 



ViX) 



sgn"(a;) (c^e'^l^l +d„e-'^l^l +o(e'"l^l)) , as \x\ ^ 00, (3.20) 



where Cn = 0, and dn is a non-zero constant. In other words, this v solution is localized. Condition 
Cn = in the asymptotics (3.20) is the condition for local bifurcations in this problem. 

Now, following the arguments laid out in Sec. 2, a scalar nonlocal bifurcation occurs when the 
?;-component of the linearised equation around the u-pulsc satisfies conditions at infinity that it 
has a purely growing component. That is, scalar nonlocal bifurcations occur when one of the vi 
and V2 solutions of the Schrodinger equation (3.8) has the following asymptotic behavior: 

i;(a;)^sgnJ'(x)(ae'^l^l+7e-'^l^l+o(e-'"l^l)) as |x| ^ 00, (3.21) 

with 7 = but some non-zero q, and j = or 1 is an integer indicating the symmetry of the v 
solution. Now, a remarkable thing happens. Because it is easy to see using the solutions of Eq. 
(3.8) obtained above, that condition 7 = in the asymptotic equation (3.21) is exactly satisfied on 
the boundary curves 

Lj = uj^^^{l3) :=n- s, 0<n-s<l, (3.22) 

where n is a non-negative integer and s(/3) was defined in Eq. (3.10). In fact, on these boundaries, 
function v has an unbounded solution 

,,=.ech..|:t^!lMii^ifMz|)! (3.23) 

fc=0 (2)^^' 

when n = 2ni is even, and an unbounded solution 

V2 = sech^x smh x ^ — TFT"! — ^ — ' (3.24) 

k=o Wkk- 
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when n = 2n2 + 1 is odd. The asymptotic behaviors of these solutions are 

v{x) sgn^{x) (^ne'^l^l + Ke^'^'^^^'^^) , |x| ^ oo, (3.25) 

where 



f 2-(|^-^s) 

4''i(i) ' n = 2ni, 



gn 



(3.26) 



\ ^ ^ ^ /no o I 1 

jw^T ^, n = 2n2 + l, 

V 2 /n2 

and ^„ is another constant which can be easily calculated. When a; < 1, the second term in 
(3.25) decays faster than e"'^'^'. Thus the coefficient 7 in the asymptotics (3.21) is zero on these 
boundaries. 

The boundary curve for scalar nonlocal bifurcations (3.22) can be written alternatively as 

/? = P^^^^ico) = i {(2n - 2u; + 1)^ - l} . (3.27) 

These boundaries are plotted in Fig. 2 as solid lines for comparison with boundaries of local bifur- 
cations (dashed lines), which are given according to (3.19) by 

P = /3^^(u;) = i {(2n + 2u; + 1)^ - l} . (3.28) 

Hence, by construction, the curves of nonlocal bifurcations simply represent the continuation of 
local bifurcation curves through co = (and mapped back up via ll> —>■ —lo, since only lo"^ appears 
in the equations). In particular, at the singular value a; = we have that ui^^^ = ■ 

Lastly, we note that the first scalar nonlocal bifurcation curve on the left of Fig. 2 needs a little 
special treatment. In fact, this solid curve is given by equation 

^ - / 7^' , -| < /5 < and a; < i 

i [1 + ^1 + 8:3], -I </3<0andu;> i. ^"^'^^^ 

In other words, the lower branch of this curve is as given by Eq. (3.22) with n = 0, but its upper 
branch is given by a different function. It can be shown that on this upper branch, the solution of 
the Schrodinger equation (3.8) also has the asymptotics (3.21) with 7 = j = 0. This curve is the 
only scalar nonlocal bifurcation boundary whose functional form is partially different from (3.22). 



3.2 Numerical results 



So we have found curves on which our proposed condition for scalar nonlocal bifurcations derived 
in Sec. 2 is satisfied. It remains to be seen what happens to the fully nonlinear equations for this 
example along such curves. In Sec. 4 below we shall consider this problem via asymptotic analysis. 
In this section we turn to numerical methods. 

First let us demonstrate further properties of the structure of solutions to the linearised problem 
(3.8) by computation of its even and odd solutions as the parameters vary. Figure 3 depicts solutions 
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of the constrained linear boundary value problems 

rX 

Vx{^) = 0, / v{x)'^dx = const. (3.30) 
Jo 

and 

rX 

v{0) = 0, / v{xfdx = const., (3.31) 
Jo 

for even and odd solutions respectively. Here X is a large positive constant. At the right-hand 
boundary point we can distinguish between solution components that decay with exponential rate 
e~^^ and those which grow with rate e^^ by considering the corresponding eigenvectors in the 
{v,v')-plaia.e. Hence we can define boundary functions 

wi =Va:{X) +ujv{X), W2 = v^{X) - u;v{X), (3.32) 

so that a zero of wi defines a solution with no exponentially growing component whereas zeros 
of W2 define solutions with no component that decays like e""^^. Hence according to the above 
definitions, w;i = can be used as a numerical test function for local bifurcations and W2 = as & 
test function for scalar nonlocal bifurcations. 

Specifically Fig. 3 depicts the results of a numerical continuation of even and odd solutions to (3.8), 
satisfying (3.30) and (3.31) respectively, for fixed a; as /? is increased from zero. It can be seen that 
an alternating sequence of zeros of functions wi and ■W2 occurs as (3 increases. At the values of each 
of the zeros we plot the mode shape v{x). Note that each successive pair of zeros corresponds to the 
function gaining an extra internal zero. The particular computation was carried out with X = 10. 
For this value it was found that the /^-values of the depicted zeros of wi and W2 correspond to 
those of the analytic formulae (3.27) and (3.28) to within 5 decimal places. Increasing X resulted 
in more accuracy, but an increase in the singularity of the boundary- value problem close to each 
zero of wi. 

For this example we have analytic formulae for the conditions defining local and scalar nonlocal 
bifurcations. Hence these computational results can be interpreted as developing numerical confi- 
dence in our method of detecting them in situations where analytic formulae do not exist (as in 
Sec. 5 below). Also, they provide geometrical insight. Thinking of the phase space (v,Vx), the 
conditions for local and scalar nonlocal bifurcations arc that the solution for large x should lie in 
one of the two eigendirections. By continuous dependence on initial condition results, then we have 
that if there are a succession of local bifurcations with increasing number of internal zeros upon 
increasing a parameter, then the solution at 'time' X must rotate in the phase plane. In so doing, 
we can not avoid having a scalar nonlocal bifurcation sandwiched between each two successive local 
bifurcations, see Fig. 4. 

Next, let us numerically investigate actual solitary wave bifurcations near the proposed scalar non- 
local bifurcation curve (3.22). First, we consider those curves with n = 1 and 2. The scalar nonlocal 
bifurcations near these curves have been numerically explored in [14]. The results are reproduced 
in Figs. 5 and 6. In each figure, solitary waves at three different locations of the parameter plane 
are shown: one is close to the local bifurcation curve (3.19) (dashed line), another one is in the 
interior, and the third is close to the theoretical scalar nonlocal bifurcation boundary (3.22). Dis- 
plays of these solitary waves are meant to show the reader how solitary waves continuously deform 
from wave and daughter wave structures as system parameters P and lo vary. As we can see, in 
both cases, scalar nonlocal bifurcations indeed occur on the theoretical curves (3.22) (see panel (c) 
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in both figures). In addition, the numerical bifurcation boundaries (circles) fall precisely on the 
theoretical curves. So for these cases, our theory is fully supported by numerics. We have also 
found similar agreement for the case n = in the nonlocal bifurcation boundary (3.22). 

However, the n = 3 case is more complicated. The bifurcation for this case is shown in Fig. 7. 
The solid line in the parameter plane (upper left panel) is the theoretical curve (3.22) for scalar 
nonlocal bifurcations. Circles are numerically detected nonlocal bifurcation boundaries. Notice that 
the numerical boundary falls onto the theoretical curve (3.22) only in the lower part. There the 
nonlocal bifurcation is indeed scalar, consistent with the geometric argument of Section 2. This can 
be confirmed in Fig. 7 (e). But in the upper and middle parts, the numerical boundary deviates 
from the theoretical curve (3.22). The reason turns out to be that, in these parts, the actual 
nonlocal bifurcation is not scalar. Indeed, an inspection of Fig. 7 (c, d) shows that the bifurcated 
solitary waves in these parts are not scalar NLS solitons pieced together. Rather, they are true 
vector solitons pieced together. Thus our analysis for scalar nonlocal bifurcations does not apply 
here. We note that near the upper part of the numerical boundary, the center of the bifurcated 
wave is a wave and daughter wave structure with n = 1 [see Eq. (3.19)], and it is flanked by two 
single-hump vector solitons on the two sides. This piecing together of different vector solitons as a 
nonlocal consequence of local bifurcation has been analytically studied in [15] before. It was shown 
there that the boundary for this type of non-scalar nonlocal bifurcation is precisely the boundary 
of local bifurcations (3.19) (here with n = 1). This is indeed the case. When the local-bifurcation 
boundary (3.19) forra = 1 is plotted as a dashed curve there, it agrees with the numerical boundary 
(circles) very well. The bifurcation in the middle part of the parameter region is also non-scalar. It 
is clear from Fig. 7 (d) that this bifurcation is somewhere in between the non-scalar bifurcation of 
Fig. 7 (c) and the scalar bifurcation of Fig. 7 (e). In fact, it is appropriate to consider this middle 
part of the bifurcation boundary as a transition between the non-scalar bifurcation in the upper 
part and the scalar bifurcation in the lower part. 

Numerical searching has revealed that there is no scalar nonlocal bifurcation observed along the 
branch corresponding to (3.22) with ra = 3 above the point at which the non-scalar bifurcations 
start. Thus the condition (3.22) can be at best a necessary condition for scalar nonlocal bifurcations. 
Why does the bifurcation deviate from scalar here? and where exactly does this deviation begin? 
These questions could not be answered by the approximate geometric argument in Sec. 2. However, 
an answer will be revealed in a matched asymptotic analysis in the next section. We will show that 
the deviation starts where the curve (3.22) with n = 3 intersects the local bifurcation boundary 
a; = 1/s of -y-pulses. 

4 Matched asymptotic theory for scalar nonlocal bifurcations in 
the CNLS equations 

To theoretically explain the scalar nonlocal bifurcation results in the previous section, an analytical 
theory will now be constructed. This theory has three objectives. The first one is to prove that 
the boundaries of scalar nonlocal bifurcations are indeed given by the condition that the solution 
of the linear Schrodinger equation (3.8) has only the purely growing component, i.e., Eq. (3.22). 
The second objective is to obtain an analytical formula for the spacing between the f-pulses and 
the central tx-pulse when the parameters are close to the boundary of scalar nonlocal bifurcations. 
The third objective is to determine when nonlocal bifurcations can deviate from scalar to non- 
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scalar. The technique we will use is similar to the tail-matching method as developed in [15] for 
the construction of multi-pulse trains, but important modifications need to be made. Throughout 
this analysis, we require a; < 1 as above. 

Suppose the ODE system (3.3) and (3.4) allows a solution where the ^-component is symmetric and 
has a dominant pulse in the center (at x = 0), while the v-component is symmetric or antisymmetric 
and has two dominant pulses on the two sides of the u-pulse (at x = ±A). Onr main assumption 
is that the f -pulses are well-separated from the central u-pulse, i.e., A ^ 1. Then we can divide 
the solution into three regions: (I) the left w-pulse region centered at x = —A; (II) the central 
it-pulse region centered at x = 0; and (III) the right u-pulse region centered at x = A. Below, 
wc will determine the solutions in each of these three regions. Note that midway between region 
II and I or III, both the u and v solutions are very small. Thus they are approximately governed 
by the linear parts of Eqs. (3.3) and (3.4), hence these solutions are linear combinations of purely 
exponentially growing and purely exponentially decaying functions to leading order. If these tail 
asymptotics from two adjacent regions can match each other, then a solitary wave can be found. 
This is the essence of the tail-matching method. 

When the u-component is symmetric or antisymmetric, the tail-matching treatment between regions 
II and III becomes the same as that between regions II and I. Thus, we will focus only on matching 
between regions I and II. 

In region I, the solution can be written as 

u = ui, V = vo + VI, (4-1) 

where 

■uo = V2ujsech[uj{x + A)], < 1, vj (4.2) 
(see Figs. 5(c) and 6(c)). In new coordinates 

e = a; + A, (4.3) 

the small ui component satisfies the linear Schrodinger equation 

ui^^ - ui + Pv^iOui = (4.4) 

to leading order. To obtain solitary waves, we demand that 

MO -^0, e ^ -DO. (4.5) 

At large positive ^ values, this uj solution must match the tails of the dominant u solution 

uo{x) = V^sechx, (4-6) 

in region II. This matching dictates that the asymptotic behavior of uj at large ^ values is 

uiiO 2V2e~^ • e^, A » ^ » 1. (4.7) 

The linear equation (4.4), together with the boundary conditions (4.5) and (4.7), completely de- 
termines the ui solution in region I. 
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Now we determine the small vj component in region I. When Eq. (4.1) is substituted into (3.4), 
and terms of order vj, vj and UjVi dropped, we find that, to leading order, vj satisfies the following 
equation 

vi^^ - u% + SvliOvi = -PujvoiO- (4.8) 

We note that it is important to retain the inhomogeneous term in Eq. (4.8), as otherwise, that 
equation with the vanishing boundary condition at negative infinity would always produce a local- 
ized solution which is impossible to match to the v solution in region II. The boundary conditions 
for solution vj are 

VI ^0, e ^ -GO, (4.9) 

and 

^ ae-^ + 7e-« - ^^^^^"'^ e^^-")^ A » ^ » 1, (4.10) 

1 — to 

where a and 7 are constants. The last term in Eq. (4.10) is contributed from the inhomogeneous 
term of Eq. (4.8). In deriving it, the asymptotic behaviors of uj and vq solutions were used [see 
Eqs. (4.2) and (4.7)]. 

Next, we determine the solutions in region II. In this region, the solutions can be written as 

u = uo + uii, v = vii, (4.11) 

where uq{x) is given in Eq. (4.6), and uij,vij <C 1 (see Figs. 5(c) and 6(c)). Here we only need to 
focus on the vu solution. This solution satisfies the equation 

viixx - '^'^vii + f3ul{x)vn = (4.12) 

to leading order. The leading asymptotic behavior of this solution at a; ^ —1 is 

1)// 7 (e""''' + <^e'"^) , -A<a;<-1, (4.13) 

where 7 and 5 are constants. If we only consider solitary waves with symmetric or antisymmetric v 
components, then the above vn solution would have the same symmetry. This symmetry condition 
would uniquely determine the coefficient 5. The constant 7 is selected by the condition that the 
tail asymptotics of vjj solution at x ^ — 1 in region II must match the v solution (4.1) at ^ S> 1 in 
region I. This matching gives 7 as = (2\/2a; + 7)e~'^^. Recall that vj ^ 1, and thus 7^1. As a 
result, to the leading order, 

7 = 2V2a;e-'^'^. (4.14) 

This matching also gives the relation 

a = 75e-'"^ = 2^/2u5e-^''^ . (4.15) 

One may wonder why the third term in the vi asymptotics (4.10) is not matched by vn asymptotics 
(4.13). In fact, there is a smaller term in the vu solution which is proportional to g^^^'^)^. This 
term arises due to the product of Uq and the e"*^^ component in the leading vu solution (4.13) [see 
Eq. (4.12)]. One can check that this term will exactly match the third term in the vi asymptotics 
(4.10). So there is no contradiction here. But this is a minor issue which is not critical to our 
analysis. 

Now we are in a position to derive a formula for the spacing A between the v-pulses and the middle 
n-pulse. This formula comes from the solvability condition for the vi equation (4.8) together with 
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the boundary conditions (4.9), (4.10) and (4.15). It is noted that Eq. (4.8) is self-adjoint, and 
it has a locahzed homogeneous solution vq^ due to the spatial translation invariance of the ODE 
equation (3.4). Calculating the integrals of products between vq^ and the two sides of Eq. (4.8) 
from — oo to y, and integrating by parts, we get 



y 

-PujVoVo^d$, = (vi^vo^ - viVoii)\^_^ ■ (4.16) 



When A » y » 1, substituting the boundary conditions (4.9), (4.10) and (4.15) into the above 
equation, we find that 



/- 



" -PujvoVo^dC -mu^e-''^^ + ^ '^ e"(i-'^)y, y ^ oo. (4.17) 

UJ 



The above equation is the leading two-term expansion for the integral on its left-hand side. When 
y approaches infinity, this integral diverges. But we can separate this divergent part from the rest 
of the integral. Notice that 



32/3w^e 

-e 



2(i-u,)y ^ n Q4pu;^e-^\^^'-^^^dC (4.18) 

J— CO 



-2A 
1-J 

Thus Eq. (4.17) can be rewritten as 

/oo r , 
ujiv^ok + 128u;3e-2^e2(i-'^)« d^ = 32Sio^e-^^^. (4.19) 
-oo '- -' 

The integral above is no longer divergent. In fact, one can use the asymptotic relations (4.5) 
and (4.7) to check that the integrand in that integral approaches zero exponentially as |x| goes to 
infinity. 

Eq. (4.19) gives a formula for spacing A when the system parameters /? and lo are specified. This 
formula can actually be made more explicit as follows. Recall that function uj is determined by 
Eq. (4.4) and boundary conditions (4.5) and (4.7). Under the notation 

%(0 = 2x/2e-^</.(0, (4.20) 

function 0(^) is then uniquely specified by the following equation and boundary conditions 

^£{-0 + /)»S(C)0 = O, (4.21) 



where vo(0 given by Eqs. (4.2) and (4.3). Under these notations, formula (4.19) simplifies as 

^-2(i-.)A ^ (4_23) 

where / is the integral 

/oo p / s n 

(/.2(t;2)^ + 16a;^e2(i-'^)« d^. (4.24) 
-oo 

Recall that the constant 5 is defined by Eqs. (4.12) and (4.13). To be more explicit, S is defined by 

ijj^^ - u'^iP + Pul{x)^ = 0, (4.25) 
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and 

V'(x) — > e-'^^ + (5e'^^, X -oo. (4.26) 

In other words, S is the coefHcient of the purely decaying component of the Schrodinger equation 
(4.25) at a; = — oo. Function tp^s boundary condition at x = oo is provided by the symmetry of 

the V component in the sohtary wave we are seeking. Since we are focusing on symmetric and 
antisymmetric v components, function would have the same symmetry. This symmetry helps to 
uniquely determine the S coefficient in the above linear problem. 

Formula (4.23) is the key result of this section. It explicitly gives the expression for the spacing 
A in solitary waves bifurcating from scalar nonlocal bifurcations. Several observations quickly 
follow from this formula. First, solitary waves from scalar nonlocal bifurcations exist only when 
parameters S and / have the same sign. Second, when 5 = 0, A goes to infinity. Thus this is a 
boundary of scalar nonlocal bifurcations. This condition is precisely the one proposed in Sec. 2. 
For the coupled NLS system, 5 = on the curves (3.22). One may notice from formula (4.23) that 
A also goes to infinity when 7 = oo. However, I = oo does not correspond to a boundary of scalar 
nonlocal bifurcations. The reason is as follows. For the coupled NLS system, 7 = oo on the local 
bifurcation boundaries of pure i;-pulses: 

u = (4.27) 
s — n 

where n is an integer and < n < s. This is because on these boundaries, the solution </> of Eq. 
(4.21) which satisfies the zero boundary condition at ^ = — oo [see Eq. (4.22)] is always localized. 
Thus in order for it to satisfy the boundary condition (4.22) at ^ = oo, (/) must be infinitely large. 
Hence I = oo. The above fact applies to the solution tt/ of Eq. (4.4) as well: on the local bifurcation 
boundary (4.27), solution uj satisfying boundary conditions (4.5) and (4.7) is infinitely large. When 
this happens, our original assumption <C 1 for scalar nonlocal bifurcations breaks down. Hence if 
there is a nonlocal bifurcation here at all, it would not be scalar: the pulses on the two wings would 
be true vector solitons. Thus / = oo does not give a boundary of scalar nonlocal bifurcations. 
Consequently, a scalar nonlocal bifurcation boundary is given entirely by the condition 6 = 0, 
which is our previous condition. An interesting and subtle issue is: what if curves S = and I = oo 
intersect? As we have discussed above, when / = oo, the nonlocal bifurcation (if there is one) 
becomes non-scalar. TIitis at the locus of I = oo and 5 = 0, the nonlocal bifurcation could turn 
from scalar to non-scalar. Then the bifurcation boundary would deviate from the scalar bifurcation 
curve 5 = at the intersection. This phenomenon could, and does, happen. In fact. Fig. 7 gives a 
good example. Let us reproduce Fig. 7's scalar bifurcation curve [i.e., (3.22) with n = 3] and the 
true bifurcation boundary in Fig. 8 below (solid line and circles). On top of it, we plot the local 
bifurcation curve (4.27) of t;-pulscs with n = (dashed line). We see that the intersection between 
these two curves is precisely where bifurcation turns from scalar to non-scalar, thus deviation 
between the numerical bifurcation boundary and the scalar bifurcation boundary starts there. This 
example and the matched asymptotic analysis above tell us that the previous condition for scalar 
nonlocal bifurcations is a necessary but not sufficient condition. When local bifurcation curves of 
f-pulses intersect with these necessary-condition curves, non-scalar bifurcations could start, thus 
scalar bifurcations on part of the necessary-condition curves will not materialize. 

The analytical expressions for 6 and / may be possible to obtain, as the linear Schrodinger equations 
(4.21) and (4.25) can be solved using hyper-geometric functions (see [7] and Sec. 3). But such 
expressions would be very complex. For practical purposes, it is preferable to determine them 
numerically. For illustration purposes, we select u = 0.6, and show the numerical values of these 
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quantities in Fig. 9 at various /3 values ranging from —1 to 7. In Fig. 9a, the integral / is shown. 
This integral is independent of the symmetry of the v solution. In Fig. 9(b, c), the coefficients 
S for anti-symmetric and symmetric v solutions are shown respectively. This figure is helpful in 
revealing on which side of the scalar nonlocal bifurcation boundary solitary waves can be expected. 
For instance, at a; = 0.6, 6 = when /3 = 0.28 and the v{x) component is anti-symmetric. On the 
left-hand side of this /? value, 6 > and J > 0, thus solitary waves with symmetric u component 
and anti-symmetric v component can be expected. But on the right-hand side of this /3 value, S < 
and I > 0, thus no such waves can be found. These predictions completely agree with the numerics 
shown in the previous section. Similar agreement is found near other scalar nonlocal bifurcation 
curves as well. 

As oj ^ oj^^^{l3), 5^0, thus A — >■ oo according to formula (4.23). Below, we derive the 
asymptotic formula for A when oj — > (P), or equivalently, p Pn^^{^) [see Eq. (3.27)]. The 

latter limit will be adopted in the following derivation as it is a little more convenient. 

When (3 — >■ f5^^^{u)), the integral / approaches a finite constant value I{lo, /3^^^), while 6 goes to 
zero. Obviously, the asymptotic formula for A crucially depends on the asymptotic formula of 6. 
We determine 6's asymptotic formula by regular perturbation methods below. Equation (4.25) can 
be rewritten as 

- a;V + Pn'^'^y'l^ = -^4^, (4.28) 
where e = (3 — Pn^^. When e = 0, this equation has an unbounded solution 

. . f —vi(x), n = 2ni, 

I ~i^2(a;), n = 2n2 + l, ^ ' 

where functions vi^2{x) and constant gn are given by Eqs. (3.23), (3.24) and (3.26). The asymptotic 
behavior of this solution is that 

V^„o(x) -^e-'"" + o(e""), x^-oo (4.30) 

[see Eqs. (3.25)]. When e is small, solution can be expanded into a regular perturbation series 

V'(x,/?,a;) = V'no(x,^™,a;) + eV„i(x,^^^^,a;) + 0(e2). (4.31) 

When this expansion is substituted into Eq. (4.25), at order e, we find that function satisfies 
the equation 

i^nlxx - W^V'nl + Pn^^ul'^nl = -^^oV'nO- (4.32) 

Recalling that the boundary conditions for ip{x) and tpno{x) are given by Eqs. (4.26) and (4.30), 
thus the boundary condition for function tpni is 

^ni{x) — ^ -e'^^, X -oo. (4.33) 

The linear operator on the left hand side of Eq. (4.32) is self-adjoint. In addition, ipno{x) is a 
homogeneous solution. Calculating the inner product between tpno and the inhomogeneous term of 
Eq. (4.32), we readily find that 



oo 

ullplodx = {ipnOxIpnl - ^n0^nlx)|!°oo " (4-34) 
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Substituting the boundary conditions (4.30) and (4.33) into the above relation, we find that for 
both i/j{x) symmetric and anti-symmetric cases, the constant S is given by the asymptotic formula 



S = ^ uUlodx + O ((/? - /3™)^) . (4.35) 

We have compared this formula with the numerical values of 6 as displayed in Fig. 9(b,c). The 
slope Uq^'^qcIx / Alo predicted by this formula is in excellent agreement with the numerical slope 
at /3 = Pn^^- When this formula is substituted into Eq. (4.23), we finally obtain the leading 
two-term asymptotic expansion for spacing function A as 

= 2i^) ^^""^^ " ^'^^''^ + lnK + 0{P- P!!^^)} , (4.36) 

where the constant K is 

K = ^n(-) = pNLBl^^^^NLB^ y_ -l^l^dx. (4.37) 



Next, we make quantitative comparisons between the spacing formulas (4.23), its leading two-term 
expansion (4.36) and numerics near the scalar nonlocal bifurcation boundaries (3.22) with n = 1 and 
2 at a; = 0.6 and various /? values (see Fig. 2). We remind the reader that a; = 1 — s is a boundary 
for anti-symmetric v components, and cj = 2 — s is a boundary for symmetric v components. At 
u = 0.6, the boundary point is at f3i^^ = 0.28 for the former case, and is at (^2^^ = 1.68 for the 
latter case. The analytical spacings from formula (4.23) and its two-term asymptotic expansion 
(4.36) are shown as dashed and dash-dotted lines in Fig. 10(a, b) for these two cases respectively 
We have also numerically determined the spacings between the ?;-pulses and the central u-pulse in 
the exact solitary waves. These numerical values are shown as solid lines in Fig. 10(a, b). We see 
that, when the separation A is large, the formula (4.23) and its asymptotic form (4.36) agree with 
the numerical values perfectly. 



5 Example 2: Saturable nonlinearity 

As a second example we take the coupled nonlinear Schrodinger system studied by Ostrovskaya 
and Kivshar [10] 

This dimensionless model arises after scaling of a model for the incoherent interaction between two 
linearly polarized optical beams in a biased photorefractive medium. Here s (0 < s < 1) is an 
effective saturation parameter, representing the photorefractive effects. This system is significant 
because it was shown in [11] that multi-humped stationary pulses may be stable solutions, a result 
that explains the experimental observations of [9]. In the limit s — 0, this system reduces to the 
Manakov equations. 
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The solitons in this system are of the form 



U{x, t) = e'%{x), V{x, t) = e'^'^^vix), (5.3) 



where u and v are real functions satisfying the following ODEs: 



Uxx-u + Uz;-. — r^- — 2\ = ^' (^■^) 
1 + s(m^ + v'') 



Vxx-'^'^v + v " ' = 0. (5.5) 

1 + s(n^ + v'') 

Looking for single-component pulses {u = 0, ot v = 0) we obtain simple planar equations, which 
can be shown by phase plane techniques to possess symmetric homoclinic orbits ztn/j, iiy'/j. Unlike 
the previous example, we know of no closed form expressions for these solutions other than at 
s = 0. So we turn straight away to numerical methods. Once again we restrict to < u; < 1 and 
look for local and scalar nonlocal bifurcations from the pulse Uh- Figure 11 shows the analogue of 
Fig. 2 for this example where the bifurcation boundaries were obtained by numerically imposing 
the boundary conditions for local bifurcation and for scalar nonlocal bifurcation on the linearised 
equation 

vxx - u^'v + = 0. (5.6) 

1 + sulix) 



Figure 12 depicts the corresponding graphs of wi{s) and W2{s) for fixed to = 0.5, where wi and W2 
were defined in (3.32). Note that we have qualitatively the same structure as in Fig. 3 in that the 
mode shapes v{x) gain increasingly many internal zeros as the parameter increases to some limit. 
Moreover, between each pair of local bifurcations of a given symmetry type there is a scalar nonlocal 
bifurcation, and vice-versa. Here it would seem that the limit s — > 1 plays the same qualitative role 
as /3 — > oo in the previous model, and we conjecture that there are infinitely many local and scalar 
nonlocal bifurcations as we approach this limit. This is bourn out by the results in Fig. 11 which 
show that the curves in the (s, u;)-plane defined by zeros of = and W2 = become increasingly 
steep as s ^ 1. It would be fair to conjecture from the numerics that all curves (apart from the 
first symmetric scalar nonlocal bifurcation) approach the point {u>,s) = (1,1). There are some 
numerical difficulties in computing the nonlocal bifurcation curves up to this point, because the 
pulse Uh becomes infinitely broad as s ^ 1, and hence the solution v{x) of the linearised problem 
which grows exponentially at infinity must be continually rescaled to avoid the solution becoming 
arbitrarily large. 

Now the presence of local bifurcations has been found before in this system; see Fig. 1 in [10] 
where local bifurcations are found via zeros of a certain integral as a function of lo for several 
different values of s. Here we have presented a simple procedure for automatically detecting these 
local bifurcations as a function of all system parameters. As far as we know these curves of local 
bifurcations are not given by closed form analytic expressions as they were in the previous example. 

The question remains whether the curves which have been putatively called scalar nonlocal bifur- 
cations really represent that for the full equations. Figure 13 depicts numerical solutions of the 
fully coupled nonlinear equations for fixed uj = 0.5 in a parameter region between the second lo- 
cal and first scalar nonlocal symmetric bifurcations. A one-parameter family of vector solitons is 
computed which are born at the local bifurcation. It is found that the solution branch terminates 
at the value s ^ 0.8055 which agrees to four decimal places with the scalar nonlocal bifurcation 



17 



boundary determined from our geometric argument of Sec. 2. Inspection of solitary waves close 
to this boundary indicates that a scalar nonlocal bifurcation indeed occurs here. Note from Fig. 
13(c) panel that the separation varies logarithmicly with the variation of the parameter away from 
the nonlocal bifurcation boundary, as in the first example [see Eq. (4.36)]. 

Finally, Fig. 14 indicates the results of a series of similar one-parameter continuations (for either 
fixed s or fixed u)) starting from the same curve of local bifurcations. Observe the agreement between 
the termination of these branches and the theoretical scalar nonlocal bifurcation boundary. Indeed, 
in each run the final orbit was found to be a pure tt-pulse flanked symmetrically by two u-pulses, 
qualitatively the same as in Fig. 13. Similar results starting from each of the local bifurcation 
boundaries computed in Fig. 11 have found that the bifurcated branches all continue smoothly up 
to the theoretical scalar nonlocal bifurcation curves. Thus, in this second example, the proposed 
condition for scalar nonlocal bifurcations and numerics show complete and global agreement. We 
note that no deviation of the sort in Fig. 7 of the first example occurs here because scalar nonlocal 
bifurcation boundaries do not intersect with local bifurcation boundaries of ?;-pulses now. In fact, 
no local bifurcations of daughter n-solutions from pure t;-pulses occur for s < 1. This is because, 
unlike in the CNLS example where the transformation cu ^ mapped back curves into the 
domain of interest (0 < /? < oo), here this transformation takes curves into the region s > 1. 

6 Linear instability of solitary waves arising from nonlocal bifur- 
cations 

An important question about solitary waves generated by scalar nonlocal bifurcations is their linear 
stability. These waves are multi-humped by construction. The common wisdom is that multihump 
solitary waves in conservative systems are linearly unstable (for instance, see [16]). However, this 
is not always the CclS©, clS stable multihump solitary waves in the saturable coupled NLS model 
(5.1), (5.2) have been reported [11]. A comprehensive stability analysis of these solitary waves is 
quite involved and lies outside the scope of the present article. However, in this section, we will 
selectively test the linear stability of a few such solitary waves as discussed above. These results are 
suggestive of the linear stability behavior of solitary waves born out of scalar nonlocal bifurcations 
as a whole. 

Our strategy for testing the linear stability is by numerically simulating the linearized equations 
of the two models (3.1, 3.2) and (5.1, 5.2) expanded around solitary waves. For both models, we 
perturb solitary waves as 

U{x,t) = e'^u{x) + e^{x,t)} , V{x,t) = e''^^^v{x) + e(l){x,t)} , (6.1) 

where {u,v) is a solitary wave, e ^ 1, and {tp,(f>) are perturbation functions. Substituting this 
perturbed solution into each model and dropping terms which are order and higher, the linearized 
equations for and (p will be obtained. Then we simulate this linearized system for a long time, 
starting with a random noise. If the solutions ip and (f) exponentially grow, then the solitary wave 
{u, v) is linearly unstable. 

First, we consider the coupled NLS model (3.1, 3.2). Its nonlocal bifurcation boundary is shown 
in Fig. 2. At {P,uj) = (0.29,0.6) which is close to the nonlocal bifurcation boundary with n = 1, 
the solitary wave born out of this nonlocal bifurcation is displayed in Fig. 15(a) [see also Fig. 5]. 
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Numerical simulation of the linearized equation around this solitary wave shows that there is a 
purely exponentially-growing unstable eigenmode. The unstable eigenfunction is displayed in Fig. 
15(b). The unstable eigenvalue is 0.083 (purely real). Thus, this solitary wave is linearly unstable. 
At {P,uj) = (1.6,0.6) which is close to the nonlocal bifurcation boundary with n = 2 (see Fig. 2), 
another solitary wave is born and displayed in Fig. 15(c). We have found that this wave is also 
linearly unstable. The unstable eigenfunction is shown in Fig. 15(d), and the unstable eigenvalue 
is 0.13 (purely real). 

Next, we consider the saturable model (5.1, 5.2). As shown in Fig. 13, at (s,u;) = (0.81,0.5), 
there is a solitary wave born out of a nonlocal bifurcation. This solitary wave is reproduced in 
Fig. 16(a). By simulating the linearized system, we have found that this solitary wave is linearly 
unstable as well. The unstable eigenfunction is shown in Fig. 16(b), and the unstable eigenvalue 
is 0.090 (purely real). 

From these selective numerical testings, we have reason to believe that solitary waves which are 
just born out of scalar nonlocal bifurcations in the two models (3.1, 3.2) and (5.1, 5.2) are linearly 
unstable in general. This belief is consistent with a previously established fact that solitary waves 
born out of vector nonlocal bifurcations are linearly unstable [16]. However, when these solitary 
waves have moved far away from the boundaries of nonlocal bifurcations, instability may disappear, 
and the solitary waves may become stable. In fact, we have found that in the saturable model (5.1, 
5.2), when solitary waves of Fig. 13(a) move to the local bifurcation boundary s = s^^ , they indeed 
become linearly stable (a similar phenomenon has been reported in [11]). But in the CNLS model 
(3.1, 3.2), solitary waves remain unstable when they move to the boundaries of local bifurcations 
[13, 12]. 



7 Discussion 



In this paper we have described by numerical computation, plausible argument and detailed asymp- 
totic analysis a new kind of bifurcation of solitary waves for coupled nonlinear Schrodinger systems. 
Our geometric argument suggests that this so-called scalar nonlocal bifurcation occurs when the 
linearization of the nonlinear system around a scalar pulse has purely growing asymptotics at in- 
finity. Our matched asymptotic analysis confirms this condition. It further reveals that when this 
scalar nonlocal bifurcation boundary intersects with certain local-bifurcation boundaries, the actual 
nonlocal bifurcation can turn from scalar to non-scalar at the intersection. All these theoretical 
results are fully supported by our numerics on two coupled NLS systems (3.3), (3.4) and (5.4), 
(5.5). Our direct numerical simulations suggest that solitary waves which are just born from scalar 
nonlocal bifurcations are linearly unstable, but they can lead to stable solitary waves by parameter 
continuation. 

It remains an open problem to fully unfold the codimension-two nonlocal bifurcation that occurs 
when the scalar nonlocal bifurcation boundary intersects the local-bifurcation boundary. In par- 
ticular it remains unknown whether a scalar non-local bifurcation continues to occur on the upper 
portion of the solid curve in Fig. 8 (upper left panel). We have found that the wave that starts 
from the local bifurcation with n = 3 does not end there, but this does not exclude the possibil- 
ity of other branches of vector solitons terminating on this possible nonlocal bifurcation boundary. 
Preliminary numerical searching has not revealed any candidate branches, but a more careful study 
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is required. 



Finally wc should mention that the ideas in this paper were originally motivated by the observation 
of what was termed a 'jump' bifurcation in [6]. There, a nonlocal bifurcation occurs where a pure-f 
pulse bifurcates at infinity from a vector soliton. The model in question is for a third-harmonic 
generation system where the symmetry is such that pure u solutions do not exist. At present it is 
not clear whether a simple criterion like the one in this paper can be applied. One could conjecture 
that the jump bifurcation occurs when linearisation around the central vector soliton has precisely 
the growing exponential asymptotics in its tail of the pure-u pulse. But this conjecture requires 
careful numerical and analytical confirmations. 
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Figure 1: Sketch figure defining (a) scalar nonlocal and (b) local bifurcations. In each case the 
top panel depicts schematically the invariant planes {u = 0} and {v = 0} and the solution of the 
linearised problem around u^- The lower two panels depict the corresponding bifurcated vector 
solitons in phase space and as graphs. 
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Figure 2: Local and scalar nonlocal bifurcation boundaries in the coupled NLS system (3.3) and 
(3.4). Scalar nonlocal bifurcation boundaries are solid lines, and are given by Eq. (3.29) [with 
n = 0] and (3.22) [with n = 1,2]. Local bifurcation boundaries are dashed lines, and are given by 
Eq. (3.19). 
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Figure 3: Curves of (a) wi{P) and (b) ■W2{P) for the linearisatfon (3.8) of the CNLS example, 
computed for u = 0.5 using and interval x G [0, 10]. Solid lines represent asymmetric solutions and 
dashed lines symmetric ones. The inserts depict solutions of the linearised equations at the first 
three zeros of wi, both symmetric and antisymmetric, which define local bifurcations; and of W2 
which define the necessary conditions for scalar nonlocal bifurcations according to our geometric 
theory. 
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(a) P = Pi 



(b) P = p2 




(C) P = |33 (d)(v(X),v,(X)) 




Figure 4: Sketch figure illustrating the large-x asymptotics of even solutions to (3.8) as the pa- 
rameter /3 varies for fixed < w < 1. The solution is depicted to the boundary value problem 
(3.30) up to x = X. Between two /3-values (a),(c) at which condition wi = for local bifurcations 
occur, there is a /3-value (b) for which W2 = 0. Panel (d) sketches the locus of boundary points 
{v{X),Vx{X) as a function of (5. 
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Figure 5: Solitary waves in the parameter region (upper left figure) bounded by the scalar nonlocal 
bifurcation boundary (3.22) (solid) and local bifurcation boundary (3.19) (dashed) with n = 1. 
Circles represent the numerical scalar nonlocal bifurcation boundary for n = 1 obtained in [14]. 
Solitary waves at stars marked by letters 'a, b, c' in the parameter region are shown with corre- 
sponding letters in the title. 
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Figure 6: Solitary waves in the parameter region bounded by the scalar nonlocal bifurcation bound- 
ary (3.22) (solid) and local bifurcation boundary (3.19) (dashed) with n = 2. Circles represent the 
numerical scalar nonlocal bifurcation boundary for n = 2 obtained in [14] . 
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Figure 7: Different types of nonlocal bifurcations appearing near the scalar nonlocal bifurcation 
curve (3.22) with n = 3 (solid Hue) in the coupled NLS system (3.3) and (3.4). Circles are nu- 
merically obtained nonlocal bifurcation boundaries. The lower part of the numerical curve lies on 
the solid line (3.22), and the bifurcation there is indeed scaler [see (e)]. The upper part of the 
numerical curve lies on the dash-dotted local bifurcation boundary (3.19) with n = 1, and the 
nonlocal bifurcation there is non-scalar [see (c)]. The nonlocal bifurcation in the middle part of the 
numerical curve is somewhere in between scaler and non-scaler [see (d)]. The dashed line in the 
upper left panel is the local bifurcation curve (3.19) with n = 3. Stars are parameter values where 
the solitary waves are shown. 
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Figure 8: Intersection between the dashed curve u = 1/s [see Eq. (3.10)] defining local bifurcations 
of daughter-u solutions from v-pulses (with n = 0) and the solid curve defining the necessary 
condition (3.22) for scalar nonlocal bifurcations (with n = 3). Circles show the numerical results 
for nonlocal bifurcations reproduced from Fig. 7. 
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Figure 9: Parameter values 6 and / at a; = 0.6 and various (3 values, (a) the integral / values; (b) 
coefficient S for anti-symmetric v components; (c) 5 for symmetric v components. 
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(b) 



w=0.6, n=2 

solid: numerical 
dashed: formula (4.23) 
dash-dotted: formula (4.36) 
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Figure 10: Comparison between the analytical formula (4.23), its leading two-term asymptotic 
expansion (4.36) and numerical values for spacing A at u; = 0.6 and various (3 values, (a) comparison 
near the scalar nonlocal bifurcation boundary (3.22) with n = 1; (b) comparison near the scalar 
nonlocal bifurcation boundary (3.22) with n = 2. 
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Figure 11: Curves of local (dashed) and scalar nonlocal (solid) bifurcations from the pure w-pulse 
Uh of (5.4), (5.5) obtained by numerically solving the linearised equation (5.6) and requiring the 
solution to have purely decaying or growing asymptotics at infinity, (a) the u-component is even; 
(b) the t^-component is odd. Only the first three curves of each type is presented. 



30 




Figure 12: Curves of (a) wi{s) and (b) W2{s) defining respectively local and scalar nonlocal bifur- 
cations according to the theory, for the model (5.4), (5.5) at w = 0.5. Dashed lines correspond 
to modes with even symmetry and solid lines to odd-symmetric modes. The labeled zeros of each 
function correspond to mode shapes for v{x) that are qualitatively similar (i.e. having the same 
structure of zeros) as the corresponding panels of Fig. 3. 
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Figure 13: (a) The t;-component of vector solitons in (5.4), (5.5) for u = 0.5 and s = 0.9075(f 



s^^), 0.9, 0.85, 0.81, 0.8055(~ s'^^^). The inset depicts the n-component of the vector sohton at 
s^^ (larger peak) and s^^^ (smaller peak). In (b) the separation between the two maxima of |f | 
is plotted as a function of s, and (c) plots the same data on a logarithmic scale where the more 
precise value s^^^ = 0.8053134 is used. 




Figure 14: A scalar nonlocal bifurcation boundary for the system (5.4), (5.5). The stars represent 
scalar nonlocal bifurcation boundaries that were obtained by solving the full nonlinear equations. 
The solid line was obtained from the condition u;2 = for the linearised equation (5.6). The dashed 
line is a similarly obtained curve of local bifurcations from which these particular vector solitons 
are born. 
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Figure 15: Linear instability of two solitary waves born out of scalar nonlocal bifurcations in the 
system (3.1, 3.2). (a, b) the solitary wave and its unstable eigenmode at {f3,Lo) = (0.29,0.6); (c, d) 
the solitary wave and its unstable eigenmode at {P,u}) = (1.6,0.6). 
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Figure 16: Linear instability of a solitary waves born out of scalar nonlocal bifurcations in the 
system (5.1, 5.2). (a) the solitary wave at {s,u>) = (0.81,0.5); (b) the unstable eigenmode. 
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